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Abstract: We present a detailed study of a filtering method based upon Dirac quasi-zero- 
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which are solutions of the euclidean classical equations of motion. On the other hand, it 
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1 Introduction 

In Quantum Field Theory the path integral is supported on distributions. Nevertheless, the 
semiclassical method aims at obtaining a reduced path integral space, labelled by smooth 
configurations. This is not in contradiction since the weights account for all the stochastic 
fluctuations around them. In the weak coupling limit this methodology works successfully 
for many theories. In the case of QCD, one might consider, for example, a dilute gas of 
instantons parametrised by its positions and sizes [1]. Classically, the distribution is Pois- 
sonian with a uniform distribution in space, a scale invariant distribution in sizes and a 
classical weight for each new instanton which is the exponential of minus its action. This is 
not the appropriate path integral measure though, since the integration over gaussian fluc- 
tuations around an instanton forces to replace the action by the free energy. This induces 
a size dependence dictated by the renormalization group. Nevertheless, such a reduced 
path-integral space, although successful in predicting a non-zero topological susceptibility, 
crucial for the solution of the U(l) problem, does not reproduce other phenomena believed 
to be present in the full theory, such as Confinement. Presumably a picture based upon 
the classical vacuum (A^ = 0) with isolated lumps is a poor description of the true vac- 
uum. Indeed, the quantum weights determined by perturbation theory already tell us that 
instantons prefer to grow and overlap with others [2, 3]. A possible modified scenario, in 
which instantons are close to each other but preserve their individuality is the instanton 
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liquid model [4]-[7]. In any case, a non-dilute picture is much harder to handle and in- 
troduces considerable complications into the semiclassical roadmap. The difficulty is not 
necessarily an impossibility. For example, the space of self-dual configurations contains 
structures with all sort of overlapping and non-diluteness, but still can be parametrised 
in terms of a few parameters known as ADHM data [8]. The dilute instanton pattern is 
included as a corner in this much richer configuration space. Other self-dual configurations 
look very different. In particular, some can be described by an ensemble of lumps of frac- 
tional topological charge 1/N. These other configurations not only lead to a non- vanishing 
topological susceptibility, but also to a non-zero value of the string tension. These and 
other attractive properties motivated a proposal for the structure of the QCD vacuum by 
our group based on fractional topological charge structures [9]. The proposal is not de- 
void of problems since a realistic description of the vacuum should allow the simultaneous 
presence of topological objects of both signs, which are not classical solutions. 

Lattice Gauge theory provides a first-principles framework to study non-perturbative 
aspects of the theory. It has been very successful in supporting QCD as the basic theory un- 
derlying strong interactions. Furthermore, it has enabled precision tests of the predictions 
of the Standard Model, a necessary ingredient for detecting evidence for any extension [10]. 
However, there are still many aspects which remain to be understood, some of which play 
a major role for candidate technicolor theories, for example. A priori, the lattice approach 
can also help in solving the questions referring to the validity of the semiclassical approach 
and the elucidation of whatever classical structures can account for the properties of the 
vacuum. Unfortunately, the general statements made at the beginning of this section imply 
that Monte Carlo generated lattice configurations contain a high level of stochastic short- 
wavelength noise. This noise completely masks the long-wavelength topological structure 
information. This justifies the necessity for finding a procedure that would enable one to 
reduce this noise and surface the underlying long-wavelength structures present in them. 

A filtering method is a procedure which takes as input a gauge field configuration with 
stochastic noise and produces a new one which preserves its long-range and topological 
structure and substantially reduces the short-range stochastic noise. Many are the methods 
that have been explored to achieve this purpose. They range from those which iteratively 
smooth the gauge field, like smearing [11]-[15], cooling [16]-[22], or the Wilson flow [23], 
to those that attempt at reproducing the topological density through a truncation of the 
spectral density of operators like the Dirac [24]- [31], or the Laplacian operators [32]. Ideally, 
any such method should produce a clean picture of the underlying long-range structures 
without introducing biases and distortions of the final data. However, all of the methods 
mentioned above have been criticized in this respect since the amount of filtering they 
introduce depends on ad-hoc parameters such as the number of cooling or smearing steps, 
or the cut-off at which the filtering operators are truncated. Still, they have been used quite 
successfully in analyzing global as well as local properties of the QCD topological structures. 
Recent reviews of the results obtained employing these methods as well as comparison 
of the performance of different filtering procedures can be found i.e. in Refs. [33]- [35]. 
The aim of the present work is to explore and develop a particular implementation of 
a filtering method, introduced in Ref. [36], which does not suffer from the ambiguities 
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mentioned above. The basic point which motivated the proposal is the idea that an optimal 
filtering method should leave smooth classical configurations invariant. Note that this is a 
non-trivial requirement, not satisfied for instance by any finite truncation of the spectral 
density of the Dirac operator in the fundamental representation. The key point exploited 
in Ref. [36] is that Supersymmetry equates the density of certain zero-modes of the Dirac 
operator in the adjoint representation with the gauge field action density [37]. However, 
in the presence of noise both observables behave in a very different fashion. On general 
grounds it is to be expected that the density of the supersymmetric zero-modes tends to 
suppress the effect of high frequency noise as compared to the action density itself, thus 
providing a filtered image of the latter. In Ref. [36] the idea was presented and a series of 
tests were carried on to show that the method works satisfactorily. 

In this paper we take this idea and develop one of its most elegant and well defined 
implementations. This can be translated to the lattice by making use of the overlap Dirac 
operator [38], which shares many nice features with its continuum counterpart. Preliminary 
results have been presented in [39]-[41]. The structure of the paper is as follows. In the 
following section we explain the definition of the filtering procedure in its continuum version 
and discuss why it is supposed to produce a filtered image of the gauge action density even 
in the presence of noise. Section 2.1 provides the details of the lattice implementation of 
the procedure in terms of the overlap Dirac operator. In section 3 we carry on a systematic 
study of how the filtering affects smooth configurations, both solutions of the classical 
equations of motion and others that are not, such as instanton - anti-instanton pairs. 
A basic goal is to quantify the distortions and biases that the method could introduce. 
The final test, presented in section 4, corresponds to the analysis of noisy configurations 
that have been obtained by applying several heat-bath Monte Carlo sweeps to a smooth 
instanton configuration. We will show that the method spectacularly succeeds in filtering 
away the ultraviolet noise. 

Having tested the method, we embark in section 5 in a possible application which is 
novel and might have important implications in hot problems such as the observed slow 
thermalization of Monte Carlo updates at high values of the coupling [42]-[44]. Starting 
with a classical configuration we apply a series of heat-bath steps at several values of the 
coupling. This introduces stochastic noise which, even for very small number of heat-bath 
steps, masks the classical configuration from which we started. The updating procedure 
also induces a motion along the space of classical configurations, changing its moduli pa- 
rameters. In particular, for instantons, which have been our main focus, the position and 
size of the instanton performs a stochastic motion. Our method can be used to study this 
motion and test if it follows the theoretical expectations. Results are quite impressive and 
might allow the numerical determination of the relative quantum weights of other families 
of classical configurations. 

The results and conclusions of the paper are summarised in the last section. 
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2 Construction of the filtering method 

As explained in the introduction, the main idea underlying the method is the relation 
between gauge field structure and zero-mode density implied by Supersymmetry. More 
precisely, we make use of the fact that given a solution of the classical equations of motion, 
the Dirac operator in the adjoint representation has a zero-mode, called the supersymmetric 
zero-mode, whose density coincides with the action density [37]. As a matter of fact, 
decomposing the mode into its two chiral components, the corresponding densities follow 
the shape of the self-dual and anti-self-dual parts of the action density. The general details 
are given in Ref. [36] and will be explained below. For the sake of clarity we will start by 
describing the method for the particular implementation which we are using in this paper. 

Given a gauge field configuration ^^(x), we can construct the corresponding Weyl 
operators D and D for the adjoint representation. The Weyl D operator is given by 
D = D^afj_, with cj^ = (I, —it) and r the Pauli matrices. The other operator D = D^^, 
is minus the adjoint of D. The Weyl matrices and ct^, satisfy the relations 

Cr^i'^u = ll^^O^a ; O'^j.CTu = rfaOoL , (2.1) 

where rja'^ and fja'^ are the 't Hooft symbols (summation over repeated indices implied). 
Using the Weyl operators we can construct the following hermitian positive operator: 

-DD = -D^D,af,a, = -r]^''aaDf,D, . (2.2) 

The spectrum of this operator is doubly degenerate. This follows from the reality of 
the covariant derivative in the adjoint representation. Given an eigenvector -0, one can 
construct an orthogonal one of the same eigenvalue as (J2''P* ■ The two bi-spinors can be 
combined into a single quaternionic eigenvector ^ = ^ ^cr^ = ("0) ^2^*)- We are identifying 
the space of quaternions with a real subspace of the space of 2 x 2 complex matrices. The 
Weyl matrices are a basis of the quaternionic space. Using this basis the operator can 
be seen as a 4 x 4 matrix whose components can be obtained by making use of Eq. (2.1). 
Having presented the basic terminology, we proceed with the explanation of the method. 

The basic object, in terms of which the filtering procedure is formulated, is O^'^^ = 
—PDDP, where P denotes the projector onto the 3 dimensional subspace generated by 
the Pauli matrices. Li components it is given by: 

0]+) = -D,D,^^'^f]f = -5jkDl - e^.krirD.D, . (2.3) 

The right hand side can be processed using the expression for the commutator of covariant 
derivatives: 

o'jV = -^3kDl + ^eijkVrF^" = -^JkDl + ieijkiE, + B^) . (2.4) 

This restriction preserves the real, hermitian, positive semi-definite character of the —DD 
operator. 
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The filtered image of the self-dual part of the action density is defined as the density 
of the eigenmode '^^^\x) of lowest eigenvalue of the operator O^^^: 

Y,0'^t^^^^\x) = X+^f-\x) . (2.5) 

k=l 

Indeed the components in colour space \I'J^^^(x) are real numbers, and the density is given 

by 

3 TV^-l 
j=l a=l 

The whole construction can be repeated for the other chirality. The main difference is 
that the ij and fj symbols are exchanged. Hence the final formula becomes 

Ojfc^ = -SjkDl - '-ei^kfirF^" = -SjkDl - ie,jk{E^ - B^) . (2.7) 

The corresponding mode ^ (x) is the filtered image of the anti-self-dual part of the gauge 
field. For the ease of notation we will be referring to the lowest eigenvectors of the O^^^ 
operators as the adjoint filtering modes or AFM modes. 

Having presented the method, we will explain the reasons why it is supposed to produce 
a filtered image of the gauge field configuration. To start with, we will focus on smooth 
configurations which are solutions of the classical equations of motion. In this case the 
filtered image should produce an undistorted image of the configuration itself. This is 
precisely the role played by Supersymmetry and the supersymmetric zero-mode. Given a 
solution of the euclidean classical equations of motion, we can decompose its field strength 
into its selfdual and anti-selfdual parts, as follows: 

Combining the equations of motion with the Bianchi identities one sees that the self-dual 
part satisfies the equation: 

r/rZ^^Xi = , (2.8) 

and a similar equation, with r] replaced by r), for the anti-selfdual part. Introducing a 
quaternionic field ^{x) = aiXi{x), one can show (using the properties of the 't Hooft rj 
symbols) that the previous equation can be written as 

D^{x) = . (2.9) 

This equation expresses that ^{x) is a zero- mode of the Dirac equation in the adjoint rep- 
resentation. This is the supersymmetric zero-mode, whose density is exactly proportional 
to that of the self-dual part of the field. Obviously, ^(x) is also a zero-mode of the positive 
hermitian —DD operator. Noticing that ^{x) is a linear combination of the Pauli matri- 
ces alone, one concludes that Xi = {^i + Bi)/2 is also a zero-mode of the operator 0'^^\ 
Now, one understands better the role played by the projection onto the three-dimensional 
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subspace leading to the construction of 0'^~^\ Indeed, the Weyl operator D might have 
many zero-modes. This is necessarily the case for classical solutions having positive definite 
topological charge (number of adjoint zero-modes > 2NQ). However, the supersymmetric 
zero- mode is indeed unique in that its quaternionic field is traceless, at all points of space 
and for all colour components. In other words it is a purely imaginary quaternionic field. 
This property, called reality condition in Ref. [36], lies at the heart of our method. We 
lack a general proof of uniqueness, but in the particular case of self-dual gauge fields, the 
ADHM construction allows us to give formal expressions for all the zero-modes. In quater- 
nionic notation they can be written as Xa<^«- The reality condition amounts to demanding 
Xo{x) = 0. This, except for some limiting cases, indeed singles out the supersymmetric 
zero-mode. Our approach in this paper has been to directly impose the reality condition 
by looking at the projected operators 0^^\ Alternatively, one could start by exploring 
the low- lying spectra of the —DD operator, imposing a posteriori the reality condition 
to select the supersymmetric zero-mode. This is the practical approach implemented in 
Ref. [36]. Although the results obtained were quite satisfactory, it requires dealing with 
the multiplicity of the —DD zero-modes and turns out to be less efficient and numerically 
more costly than the present proposal. 

What do we know about the excited states of the spectrum of the O*-"*"^ operator. In 
general, for a finite volume, we expect a gap to develop in the spectrum. The gap goes 
to zero for large volumes, but the corresponding limiting states become non-normalizable. 
The operator has a single point discrete spectrum and a gapless continuum spectrum. 
Furthermore, the eigenmodes would look rather structureless and easily distinguishable 
from the supersymmetric zero-mode. We will see how these expectations are realized in 
our numerical results. 

Certainly a possible analytical testing ground is that of a single SU(2) instanton. Here 
the only scale is clearly p, so that we might use it as a unit of lengths and energies. The 
operator can be rewritten using the formula of the covariant Laplacian for an instanton 
in SU(2) [2] 




(2.10) 

In the previous formula: 

and i\ = —(x^dy — Xyd^Y /S, where summation over repeated indices is understood. 
In addition T and T' are the generators of SU(2) in the adjoint representation acting 
respectively on the colour or on the spinorial indices. 

If we restrict to purely radial eigenvectors and split upon eigenstates of the direct 
product representation T + T', with isospins i = 0, 1, 2, we get: 

= -(^\-^-^+ ^ + ^ t{t + l)-6 

\drj r 9r (r2 + 1) ^ (r2 + 1)2 2 ' ^' ' 
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We might write the rotationaUy invariant eigenvectors as: 

(r2 + 1) 



* = 7^^ • (2-12) 



The eigenvalue equation then reduces to 

-xG"(x) + 2|^G'(x) + |^G = ^G(x) . (2.13) 

Obviously for t = 0, G{x) = 1 is a solution of vanishing eigenvalue. This is the expected 
supersymmetric zero-mode. It is easy to see that there are in addition 5 non-normalizable 
zero-modes for t = 2 and G{x) = (r^ -f- 1)^. 

Up to now, we have just analysed the behaviour of the filtering procedure for solutions 
of the classical equations of motion, for which the filtering mode is just the supersymmetric 
zero-mode. For other smooth configurations the lowest eigenvalue will not be zero. Indeed, 
if we take the self-dual part of the gauge field Xii^) = {Ei{x) + i?j(x))/2 and apply the 
operator to it, we get 

0^^\k{x) = -\^^''D,f , (2.14) 

which only vanishes if the current, = D^F^'^, vanishes. From here we get an upper 
bound on the minimal eigenvalue of 0^~^^: 

+ ^ fdxTT{f{x)f{x)) 

- AjdxTr{xkix)xk{x)) ' ^ ' ^ 

If the smooth configuration is made of local structures which are approximate solutions 
of the classical equations of motion, the current would be very small in those regions and 
the filtered mode ^'^^^(x) will be close to the self-dual action density. The limit to this 
situation would be the region in which level crossing appears. We will explore this situation 
in our numerical analysis, to be presented later. 

In order to understand the behaviour of the O^^^ operators for non-classical configu- 
rations, one could study small perturbations about the classical solutions: + e6A^. To 
derive the appropriate formula it is better to regard O^"*"^ as the projected version of —DD. 
The eigenvalue equation, to leading order in e, takes the form 

- PD5D^ - PDD5^ = JA^- . (2.16) 

Projecting onto ^ one realizes that 5X vanishes (hence it is order e^). The equation for the 
eigenvector is 

PDD5^ = -PD6D^ . (2.17) 
The right-hand side can be evaluated to 

iauapaj[D^5Ap, ^j] + 2iaj[5Ap, Dp^j] . (2.18) 

Now, without loss of generality we might impose the background gauge condition on 5Ap. 
With this we can rewrite the right-hand side as 

iakaj [5Ek + SBk, ^j] + 2i<jj [8Ap, Dp^^] . (2.19) 
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After several operations one arrives at the equation 



- Dl6^j + ie,Jk[E^ + B^, 6^k] = -^e^jk[SEi + 6Bi, ^k] - 2i[6A^,D^^,] , (2.20) 

where = {Ej + Bj)/2. 

To see the filtering properties of the method, we might consider a deformation 6A^ 
involving high frequencies. For that purpose we make 6A^ proportional to e*^^ for high 
values of p. The corresponding field strengths 5Ei, 6Bi will then be proportional to pe*^^. 
Thus, in the corresponding action density the noise terms are amplified with a factor p^. 
Let us now look at the expected behaviour of the filtered mode S^k- The right-hand 
side of Eq. (2.20) will also contain terms proportional to pe*^^. To solve the equation 
we should decompose this term into eigenstates of the O^"^) operator. The decomposition 
should be dominated by states of high eigenvalue, which are essentially given by plane 
waves e*''^. Since, the field-strength of the zero-mode is smooth, the typical wavelengths 
involved are g ~ p. The same reasoning tells us that the corresponding eigenvalue of the 
0(+) 

operator is q^. Applying, these considerations to Eq. (2.20) we conclude that 6"^ goes 
like p/p'^ = l/p. The new filtered mode density has a 1/p suppression of the added noise. 
This is precisely the desired effect, which makes our method a filtering procedure for the 
gauge field configurations. 

2.1 Lattice implementation 

There are two important ingredients involved in the construction of the O^^^ operators: 
chirality and the reality of the covariant derivative. In implementing the method on the 
lattice it is convenient to preserve these properties. This can be achieved by using the 
overlap operator [38] as lattice Dirac operator. This will be explained below. Some results 
obtained with this procedure were presented in Ref. [40]. 
The lattice version of the O^^^ operator is given by 

O^^^ =PP±ij5DifP±P , (2.21) 

where P± = |(lib75) and P is the projection onto the space of purely imaginary quaternions 
as before. The overlap Dirac operator is given by the formula: 

Di = —{I + j,e) , (2.22) 
a 



where e is a hermitian, unitary matrix given by: 

H 



(2.23) 



and H = 75 (—Af -|- aDy^^). The matrix D^^^ is the standard massless Wilson-Dirac 
operator in the adjoint representation given by: 



aD^jjin, m) = ^ ((I - 7^)C/^(n)5mn+^ + (I + 7/,)[/^^(m)5^„_^ - 2I5„, 



(2.24) 
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where n and m label points in the 4-dimensional lattice of spacing a, and fi is the unit vector 
in the n direction. The gauge field links in the adjoint representation can be obtained in 
terms of links in the fundamental as follows: 

{U^)^^ = 2Tr(ulTaU,n') , (2.25) 

where Ta are the generators of SU(N) in the fundamental representation with the standard 
normalization TT:{TaTi,) = ^Sab- Notice that the links in the adjoint representation are real 
orthogonal matrices. 

The overlap Dirac operator is known to satisfy 75-hermiticity and the Ginsparg- Wilson 
relation 

b5^Di} = ^Di^5Di . (2.26) 

These two conditions imply that (75^?^,)^ commutes with 75. The positive and negative 
chirality components are the lattice counterpart of —DD and —DD. 

Finally, we should verify that the projection onto the purely imaginary quaternions 
operates in the same way as in the continuum. The reality of the gauge field links in the 
adjoint representation implies the double-degeneracy of the eigenvalues of H. This follows 
from the equation 

H* = C^^HC , (2.27) 

where C is the charge conjugation matrix satisfying 7* = C~^^^C . Eq. (2.27) also holds 
for the overlap Dirac operator and for (75-Doi))^- In the Weyl representation of the Dirac 
matrices 

the matrix C can be chosen as 757072- Thus, following the same procedure as in the con- 
tinuum, the chiral components of (75-Doi))^ can be seen as a matrix acting on quaternions. 
Restriction to the subspace of purely imaginary quaternions leads to the O^^'^ matrices 
defined in Eq. (2.21). 

In our lattice implementation of the overlap operator, we have used the rational ex- 
pansion for the e matrix proposed in [45, 46]: 

„ n 

6 = hm - V , (2.29) 

n-,00 n ^ + as 

s=l 

where = cos~^(^(s — ^)) and as = tan^(^(s — |)). The action of the (if^ + Us)"-'^ 
matrices over a vector x is computed by using the method of multiple shifts described in 
[47]. As a compromise between performance and stability we have taken n = 60 and the 
tolerance parameter that controls the conjugate gradient accuracy is taken to 10~^. In 
order to minimize the storage needs, we use the two-pass algorithm proposed by Neuberger 
[48]. The pseudo-code necessary for the implementation is well described in Ref. [48]. 

In order to compute the low-lying spectrum of the O^^^ operators we have used the 
conjugate gradient algorithm [49]. The eigenvectors are computed with an accuracy such 
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that: 



,(±) 



L 



< 00 



(2.30) 



with uj = 10^®M^/a^. This guarantees no appreciable errors in the eigenvector densities. 
The speed at which convergence is achieved depends very much on the gauge field con- 
figuration and on the initial vector for the conjugate gradient algorithm. Typically, for 
the Q = I cases analysed in the paper, the computer time needed to determine the lowest 
eigenvector grows as A'^^, with a proportionality constant that increases as /3 decreases. 
In some cases we have seen that a suitable choice of the initial vector for the conjugate 
gradient can reduce this time by a factor of 5. 

To check that our results are not biased by the choice of the e matrix expansion, we 
have repeated the computation of the low-lying spectrum on some of the configurations 
with the Lanczos algorithm proposed in [50]. The results obtained by both methods are 
compatible within errors. 

Finally, we comment on the choice of the Wilson-Dirac mass term M. Most of our 
results were obtained for M = 1.4. This value is a compromise which allows the filtering 
method to detect fairly small instantons without sacrificing performance [39, 41]. Never- 
theless, we have done several tests to check the dependence of our results upon the value 
of M. In all cases, including noisy and smooth gauge field configurations, no significant 
changes are observed on the density of the AFM. This contrasts with the results in Ref. [31] 
obtained using the modes of the Dirac operator in the fundamental representation as filter- 
ing modes. In that case, the claim is that the value of M monitors the amount of filtering 
obtained for the configuration. 

3 Testing the method with smooth configurations 

In this section we will test the numerical implementation of the filtering method on a set 
of smooth lattice configurations. We will focus on two generic situations. First, we will 
look at lattice configurations that correspond, in the continuum limit, to solutions of the 
classical equations of motion. This includes discretized Q = 1 instantons and calorons 
[51, 52] as well as higher topological charge objects. The filtering method should work well 
in this case, but the analysis will allow to quantify discretization and finite volume effects. 
Second, we will analyze smooth configurations which are not solutions of the continuum 
equations of motion. We will focus in particular on the case of instanton - anti-instanton 
pairs, testing the ability of the filtering method to detect structures in both chiral sectors. 
In all cases we will analyze the low lying spectrum of the oj^^^ operators, paying particular 
attention to the occurrence of level-crossing that could disturb the identification of the 
adjoint filtering mode. 

3.1 Classical configurations 

The filtering operator is designed to select the supersymmetric zero-mode if applied to 
classical backgrounds. This is guaranteed in the continuum but, for the method to work, it 
is essential that this property is preserved when both the O^^^ operators and the classical 
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configurations are discretized. Lattice artefacts can be considered as a particular case of 
the deformations discussed in section 2. We expect them to induce an O(a^) lifting of the 
AFM eigenvalue and a small modification of its density through mixing with the excited 
states. Although this contamination goes away in the continuum limit, it is important to 
quantify the distortions it may introduce in the procedure. 

In what follows we will analyze several Q = 1 and Q = 2, SU(2) configurations dis- 
cretized on the lattice. We will show that the method spectacularly succeeds in identifying 
classical lumps down to sizes of the order of two lattice spacings. Results will be pre- 
sented for SU(2) instantons but we have also checked that the filtering procedure works 
well for other classical solutions as finite temperature calorons [51, 52], or fractional charge 
instantons [53, 54]. 

3.1.1 Single instanton configurations 

For the purpose of the test we have generated a set of Q = 1, SU(2) lattice instanton config- 
urations for several lattice volumes {Nf) and two different sets of boundary conditions. In 
addition to a periodic torus (p.b.c.) we have also considered a torus with twisted boundary 
conditions (t.b.c.) [55], corresponding to m = (0,0,0) and k = (1,1,1). The reasons to 
employ twisted boundary conditions are twofold. To start with, the gap in the spectrum 
of the filtering operator is strongly dependent on the boundary conditions. It is easy to 
understand why this is the case by analysing the spectrum on a trivial background. For 
t.b.c. the set of classical vacua is discrete and the spectrum has a gap of size vr^/L^, for 
= and finite box size L = Nga. Periodic boxes, on the contrary, have a continuum of 
classical vacua [56, 57] which support a continuum of constant density adjoint zero-modes. 
On the background of instanton configurations, these vacuum zero-modes become quasi- 
zero- modes for p « L and can give rise to unwanted level-crossing in the spectrum. This 
situation is avoided in the twisted case, where the tt^/L^ gap in the spectrum pushes any 
level-crossing to very large values of L. There is a second reason to select twisted instead of 
periodic boundary conditions for the analysis of Q = 1 instantons. On a finite periodic box 
there is a restriction over the allowed topological charge of classical solutions, and Q = 1 is 
not feasible [58]. Although this restriction disappears in the thermodynamic limit, it gives 
rise to finite size effects that can shift considerably away from zero the lowest eigenvalue 
of the filtering operator. This limitation disappears, however, as soon as one imposes a 
non-zero twist. 

The results we will present correspond to instanton configurations of various sizes 
generated by cooling, an iterative update process that monotonically decreases the lattice 
action [16]. Since the lattice breaks the scale invariance, the fate of instantons under cooling 
depends crucially on the choice of lattice action employed in the cooling algorithm [20]. One 
can take advantage of this and use cooling to modify at will the instanton size parameter 
p. In this paper we have made use of the e-cooling method proposed in Ref. [20]. It 
employs a variant of the Wilson action containing 1x1 and 2x2 plaquettes with relative 
coefficients tuned to control the 0{a'^) scale violations of the lattice action. In particular, 
e = is the choice that eliminates these corrections. Figure 1 displays the dependence of 
the integrated e = instanton action versus the instanton size for several lattice volumes. 
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Figure 1. We display the dependence of the mstanton action (in units of Stt^) versus the instanton 
size (in lattice units) for several values of the lattice volume: Ng = 14, 16, and 20. (Left) Comparison 
of twisted and periodic boundary conditions; (Right) A blow up of the plot for twisted boundary 
conditions (note the change of scale). The action corresponds to the discretized e = lattice action. 



For t.b.c. the deviation from the continuum instanton action amounts to less than one part 
per mille for instantons of size larger than p = 2.5 a. In that range we expect very small 
distortions in the spectrum of the O^^^ operator. The situation is not as good for periodic 
boundary conditions since, as we have just discussed, there are no exact classical solutions 
with Q = 1 on a finite periodic box. For this reason the deviation from the continuum 
action is generically much larger than with t.b.c. The results that will be presented in this 
section correspond to several representative twisted and periodic configurations selected 
among the ones displayed in these plots. 

For all the configurations analyzed we have computed the four lowest eigenvalues in 
each chiral sector. The sector of negative chirality is structureless. For positive chirality, 
as expected, one of the eigenvalues is distinctively smaller than the others. Its eigenvector 
provides the adjoint filtering mode (AFM) which should correspond to the supersymmetric 
zero-mode whose density reproduces the instanton action density. A comparison between 
the density of the AFM and the self-dual part of the action density is displayed in figure 2 
for two instantons of sizes p = 3a and p = 6.2a, with Ns = 20 and twisted boundary 
conditions. The qualitative agreement is impressive. For a more quantitative comparison 
we give in table 1 the values of the instanton parameters obtained from the densities of 
the AFM (Xafm, Pafm) and the self-dual part of the gauge field strength {X, p). They have 
been extracted by fitting the 1-dimensional instanton profiles to the data, following the 
procedure described in the appendix. Differences in the extracted instanton locations are 
negligible. As for the sizes, the differences remain much smaller than those associated to 
different fitting procedures - see figure 15 in the appendix. 

We proceed next to present a detailed analysis of the effect of the numerical implemen- 
tation on the spectrum. It is important for the success of the method that no level-crossing 
takes place between the lowest eigenvalue and the excited states since this would induce a 
misidentification of the supersymmetric zero-mode. We will discuss first the case of twisted 
boundary conditions which from this point of view provides a much cleaner situation. As 
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Figure 2. Left: the self-dual part of the action density (in units of Stt^), and Right: The density 
of the lowest eigenvector of the O^^'' operator, both integrated over the two orthogonal directions. 
The top and bottom plots correspond respectively to Q = 1 instanton configurations of sizes p ~ 3 a 
and p = 6.2 a. All the data are for twisted boundary conditions and Ng ~ 20. 
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Table 1. For two different size Q = 1 instantons, on a 20"* twisted lattice, we compare the 
instanton parameters determined from the 1-d profiles of the AFM density and the self-dual part 
of the action density. The corresponding size parameters are given by: Pafm and p. The difference 
in the position of the instanton center coming from both determinations is given by: AX^. 

discussed before, the spectrum of the O^^^ operator is gapless in the infinite volume limit 
with a continuum spectrum that corresponds to non-normalizable states. For finite volume 
however there is a finite gap of order vr^/L^ and no level-crossing is expected. We will an- 
alyze in what follows how lattice artifacts modify this expectation, and we will see that in 
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Figure 3. For a set of Q = 1 twisted instaiiton configurations of size p, we display: Left: the 
eigenvalue of the AFM, as a function of p/a. The line corresponds to the fit in Eq. (3.1). Right: the 
first excited eigenvalue of the o'\^^ operator, as a function of p/L. The line is the fit in Eq. (3.2). 



the worst case situation level-crossing will only occur for very large lattice sizes Ng > 500. 
Figure 3 displays, for twisted boundary conditions and two different lattice sizes, Ng = 16 
and 20, the lowest (left) and first excited (right) eigenvalues of the oj^^^ operator. Let us 
first focus on the analysis of the lattice artifacts on the AFM eigenvalue (Aq). In the plot 
it is displayed as a function of the instanton size in lattice units: p/a. The line drawn in 
the plot corresponds to a fit of the form: 



p^Xo = a{^)\f^{^y , (3.1) 



.p/ \p 

with coefficients a = 2.0(2)10"^, and /3 = 44(2) 10~^. It describes well the data for p > 3a 
and shows that the eigenvalue of the filtering mode approaches zero in the continuum limit, 
as expected. Results are given for instantons of size larger than p = 2a to avoid problems 
in the convergence of the conjugate gradient algorithm. 

With respect to the first excited eigenvalue (Ai) a comparison between the Ng = 16 
and the A^^ = 20 results shows that it scales as 1/L^ . The data can be fitted to 



where b = 3.8(1), and c = 0.16(1). From the fits to the two eigenvalues we can extract the 
ratio Ao/Ai and estimate when it will become of order 1. We find that for instantons of 
size p = 2a (3a) the level-crossing will take place for extremely large values of the lattice 
volume, Ns ~ 500 {Ng ~ 3000). 

The situation is much worse for periodic boundary conditions since in that case the 
gap in the spectrum goes to zero with decreasing p/L. The results are presented in figure 4, 
where we display the two lowest eigenvalues of the filtering operator. The AFM eigenvalue, 
displayed on the left plot, matches the twisted result, Eq. (3.1), for small p/L. It becomes 
however strongly affected by finite size effects for p > L/5, as a consequence of the absence 
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Figure 4. For a set of Q = 1 periodic instanton configurations of size p, we display: Left: the 
eigenvalue of the AFM, as a function of p/a. The line corresponds to the fit in Eq. (3.1). For 
comparison we also include the data for twisted boundary conditions and Ns = 20. Right: the first 
excited eigenvalue of the O^^'' operator, as a function of p/ L. The line is the fit in Ec^. (3.3). 



of exact Q = 1 solutions in the periodic case. The first excited eigenvalue, displayed on 
the right plot, goes to zero as /L'^- The fit in the figure corresponds to: 

with a = 1.2(1). Combining this fit with the one describing the lattice artifacts Eq. (3.1), 
we can again estimate when the ratio between the two eigenvalues becomes of order 1. In 
this case it happens for rather small values of N^: for p = 2a (3a) the level-crossing will 
take place at Ns ~ 20 {N^ ~ 35). 

3.1.2 Two instanton configurations 

In this subsection we will study other classical solutions. As in the previous case, the filter- 
ing method is bound to work in the continuum, so that our goal is mostly that of exploring 
possible problems associated to the discretization. This is a necessary step previous to 
any application to Monte Carlo generated configurations. Our study has focused in one 
situation which might potentially cause problems. This occurs whenever the gauge field 
configuration contains widely separated structures. It is to be expected that in this case 
the gap narrows down, since every isolated structure should give rise to a quasi-zero-mode 
of the O^j^^ operator. There is still a unique zero-mode associated to the super symmetric 
zero-mode, but lattice corrections might induce unwanted mixing, whose effects we want 
to explore. 

For definiteness we have concentrated our study in the case of two instantons (I-I) at 
variable separations. For the analysis we have prepared a set oi Q = 2 SU(2) instanton 
configurations, by applying the following engineering techniques. Typically, we start with 
a couple of instantons on a lattice of size 16^. The configurations are glued along the time 
direction, and subsequently cooled until the total action drops to a value close to S = IGvr^, 
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Figure 5. For a set of Q = 2 instanton configurations, on a 16^ x 32 periodic lattice, we 
display the three lowest eigenvalues of the O^^"* operator as a function of the inverse I-I over- 
lap: d/p (distance over average size, p, of both instantons). The line corresponds to the fit 
Aia2 = 0.036(5) exp{-l. 7(1) d/p}. 

signalling that the configuration is approximately a solution of the classical equations of 
motion. The gluing can be done at variable separation among the two instantons, and this 
is not substantially altered by the cooling process. Finally, we end up with a collection of 
Q = 2 configurations on a 16^ x 32 lattice and with periodic boundary conditions. Notice 
that, in contrast to the Q = 1 case, in this case there are indeed classical solutions in the 
continuum living in a box without twist. 

We studied the low-lying eigenstates of the O^^^ operator for the afore-mentioned 
set of configurations. Our results for the spectrum, displayed in figure 5, agree nicely 
with our expectations. There is a unique zero- mode (up to lattice corrections) for all I-I 
separations. The first excited eigenvalue, however, depends strongly upon separation d, 
dropping to zero as this separation gets large. This dependence can be fitted to Xia^ = 
0.036(5) exp{ — 1.7(1) d/p}, where p is the average size of the two instantons. The second 
excited eigenvalue, on the contrary, does not depend significantly upon the separation d 
and is presumably mostly determined by box-size and boundary conditions. 

Our interpretation of the results is substantiated by the structure of the two lowest 
eigenvector densities, displayed in figure 6 for several values of the separation d. The 1-d 
profiles represent the integral of the densities over the spatial directions. For comparison, 
the spatial integral of the self-dual part of the action density is also displayed. Two regimes 
are clearly observed. For an I-I distance oi d = 1.6p, the situation is similar to the one 
instanton case: The ground state density matches the self-dual action density, and the first 
excited state corresponds to a nearly constant mode. For larger separations a level crossing 
takes place in the excited spectrum and both densities display the distinctive two instanton 
profile of the action density. 
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Figure 6. For a set of Q = 2 instanton configurations, on a 16'^ x 32 periodic lattice, we display 
from left to right the one-dimensional integrated profiles of: the self-dual part of the action density; 
the densities of the lowest and first excited eigenvectors of the O^^'' operator. All the densities are 
normalized in the same way. The plots from top to bottom correspond to I-I separation of c? = 1.58p, 
2.03p, 2.72p, 3.22p, and 3.75p. The parameters for the lines in the plot have been obtained from 
the fits to the 1-d continuum instanton profiles described in the appendix. 

Let us examine these results in the light of our expectations. In the limit of infinite 
separation each instanton (n = 1,2) has an associated supersymmetric quasi-zero-mode 
Xi"^ = E^"\ given by its corresponding electric field. As the instantons approach each 
other, the ground state of 0^^\ which is an exact zero-mode, is given by the symmetric 
combination '^f'^ oc E^^^ + Ef'\ The antisymmetric combination should correspond to 
the first excited state. Both densities coincide with the self-dual action density (except 
for minor corrections in the region where the instantons overlap). On close observation, 
however, we notice that the density profiles on figure 6 show slightly different weights for 
each of the two peaks in both the ground state and the excited state. This is presumably 
an effect of the discretization. Lattice artefacts induce a mixing of both states, which can 
be estimated by diagonalizing the O^^^ operator in the subspace spanned by '^^^^ and 
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Table 2. We give the parameters of the set of Q = 2 instanton configurations analyzed in the 
paper. All the data correspond to a 16'^ x 32 lattice with p.b.c. The instanton sizes, in lattice 
units, arc extracted from a fit to the, spatially integrated, 1-d profiles of the AFM (jOafm), the first 
excited state (pai), and the self-dual part of the action density (p). The separation between the 
two instantons is given by d. By p we denote the average size of the two instantons. The mixing 
between the two lowest modes of the 0'+) operator is given by tan0. 

Writing the perturbed states as: 

§(0) = cos6'^'(°) -sin(9*(^) , 
^'(^^ = sin6l^'(°) +cos6'*(^) , 

one obtains 

tan6' = °^ , 

AA + ^AA2 + 4(505t^)2 

where AA is the gap between the unperturbed states. The lattice artefacts enter through 
the expectation value: 50^+^ = {m^^'^\50^j^\^^^^). We can estimate the size of tan 6 by 
replacing in the denominator of this expression the actual gap, extracted from the data in 
figure 5, and for the numerator the typical size of the O^^^ lowest eigenvalue on a Q = 1 
instanton background. For our I-I configurations we estimate that the maximal mixing 
is attained for the largest separation d = 3.75p, being of order tan^ = 0.05. A direct 
determination based on the relative heights of both peaks gives a value of tan^ = 0.04, in 
remarkable agreement with our crude estimate. 

Although the possible alteration of the peak heights, for the case of sufficiently isolated 
structures and corresponding small gaps, is a matter of concern, we want to emphasize that 
it is not an insurmountable difficulty, once the origin of the problem is pinned down. Indeed, 
our determination of instanton parameters is largely insensitive to this normalization. This 
is clearly shown in table 2, where we collect the size parameters p^") of both instantons 
extracted from fitting the action density and the two lowest eigenmode densities. As 
described in the appendix, the fitting has been done using the shape around the peaks 
of the 1-d profiles obtained by integrating the densities in the spatial directions. The 
continuous lines in figure 6 correspond to a superposition of the Id-instanton profiles given 
by Eq. (A. 3) in the appendix, with parameters, including normalization, determined from 
the fits. 

To summarize, the filtering procedure succeeds in reproducing the Q = 2 action density. 
We have studied a distortion in relative normalization of the profiles induced by lattice 
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Figure 7. For a set of I-A configurations, on a 12^ x 24 periodic lattice, we display the two 
lowest eigenvalues of the O^^'' operator as a function of the inverse I-A overlap, given by 1/0/a = 
2d/{pi + pa)- The line is a fit to Eq. (3.4). 



artifacts, which is present when the instanton separation is large enough. There might be 
several ways to cope with it, once we know its presence and origin. In particular, possible 
improvements of the O^^^ operator, which reduce the size of lattice artefacts, will also 
reduce this effect. 

3.2 Smooth non-classical configurations 

Having tested successfully the performance of the filtering method on classical configura- 
tions, we now turn to the analysis of smooth backgrounds that are not solutions of the 
classical equations of motion. The prototype example is an instanton-anti-instanton (I-A) 
pair which approaches a classical solution only in the limit of infinite I-A separation. Al- 
though in this case we do not expect exact zero-modes of the O^^^ operators, we will show 
that one eigenvalue in each chiral sector remains in general significantly smaller than the 
others. This is so even for configurations in which the instanton and the anti-instanton 
overlap considerably and allows to unambiguously select the adjoint filtering modes. As 
we will see, the densities of the AFMs correctly reproduce the self-dual and anti-selfdual 
parts of the action density. 

For the analysis, we have generated several SU(2) I-A congurations with different I- 
A separation on a 12'^ x 24 lattice with periodic boundary conditions. The set has been 
obtained by letting a well separated I-A pair annihilate under e = cooling. The initial 
configuration is engineered by gluing along the time direction an instanton with the anti- 
instanton obtained from it by time reflection. In the process of destruction of the pair the 
action monotonically decreases from values close to IGvr^ to less than Svr^, a point after 
which one can no longer identify the I-A lumps in the gauge action density. 
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Figure 8. For a set of I-A configurations, on a periodic 12'^ x 24 lattice, wc display: Left: the 
profiles of the self-dual (red) and anti-self-dual (blue) parts of the action density integrated over 
the three orthogonal coordinates, and Right: the same for the densities of the AFM modes on the 
positive (red) and negative (blue) chirality sectors. All the densities are normalized in the same 
way. From top to bottom the plots correspond respectively to inverse I-A overlap: 1/Oia = 3.35, 
2.42, 1.87, and 1.59. The lines in the plot correspond to fits to the continuum 1-d instanton profiles. 

Let us present first the analysis of the O^^^ spectra. Figure 7 displays the two low- 
est eigenvalues of the O^^^ operator as a function of the inverse I-A overlap: 1/0/a = 
2d/{pi + pa)- a similar picture is obtained in the negative chirality sector. The lowest 
eigenvalue increases with I-A overlap but a clear gap is maintained as far as Oja < 1/1.6. 
This allows to clearly identify the AFM. We will discuss the closing of the gap later on but 
let us concentrate now on the AFM densities. They are represented in figure 8 for several 
characteristic stages in the process of I-A annihilation corresponding to I/Oja = 3.35, 
2.42, 1.87, and 1.59. We present on the left the 1-dimensional profiles of the self-dual and 
anti-selfdual parts of the action density compared with the AFM densities displayed on the 
right. The agreement is in all cases very good, even for the lowest plot that corresponds to 
the maximal I-A overlap in figure 7, above which level crossing takes place. We extract the 
(anti) instanton parameters by fitting these 1-d profiles following the procedure described 
in the appendix. The results are presented in table 3, and the corresponding fits repre- 
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Table 3. Wc give the parameters of the set of instanton-anti-instanton configurations analyzed in 
the paper. All the data correspond to a 12'^ x 24 lattice with p.b.c. The instanton sizes, in lattice 
units, are computed from the 1-d density profiles of the positive (Pafm) negative (Pafm) chirality 
AFM modes, and the self-dual (p^) and anti-self-dual (p^) parts of the action density. The I-A 
overlap is given by Oia- 

sented by the continuum lines on the plot. Even for the extreme case with largest overlap, 
the difference between the size parameters extracted for the AFM modes and the gauge 
densities amounts only to 10%, decreasing down to less than 1% for Oja = 1/3.35. 

One final remark is that, unfortunately, the level-crossing is controlled not only by 
the growth of the lowest eigenvalue Aq with the I-A overlap, but also by the finite volume 
dependence of the first excited eigenvalue Ai. We can estimate in general when it takes 
place. The lowest eigenvalue can be fitted to 

fi2Ao = 21(l)exp(-^) , (3.4) 

where the fit corresponds to the continuum line in figure 7. The first excited eigenvalue 
is expected to behave as 1/L^ for twisted, and as p^/L^ for periodic boundary conditions 
(see respectively Eqs. (3.2) and (3.3)). By requiring Aq < Ai we obtain: 

for p.b.c. , (3.5) 

for t.b.c. . (3.6) 

For example, the identification of a pair with d = 1.6p requires L < 2d for periodic 
boundary conditions, and L < 5d for the twisted case. Once again the use of twisted 
boundary conditions pushes away the region in which level crossing takes place. 

4 Testing the method with noisy configurations 

We now come to the crucial test on the filtering method. Our ultimate goal is to apply 
the method to Monte-Carlo ensembles and extract from it a clear picture of its underlying 
long-range topological content. The only way in which one could trust the results is if 
one has previously tested the method in a controlled situation. It is not just a question 
of recovering the long-range information stored in the configuration, but also of making 
sure that the filtering procedure does not introduce important distortions in the main 
parameters describing the underlying topological structure. 
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Figure 9. We display from left to right: the self-dual part of the action density (in units of Stt^) 
of a smooth Q ~ 1 twisted instanton configuration, the same quantity after applying Nh heat bath 
sweeps to it, and the density of the AFM on the heated configuration. All quantities are integrated 
in two dimensions. The top and bottom plots correspond respectively to: — 14, = 10, and 
P = 3] N, = 16, Nh ^ 60, and (3 = 8. 

Our proposed strategy is the following. Starting from a smooth configuration, we begin 
to apply Monte Carlo updating steps to it. We limit ourselves to large values of /? and a 
small number of Monte-Carlo heat-bath sweeps to guarantee that the initial configuration 
is not destroyed in the process, and that only statistical noise is added to it. In this and the 
following section we will test the ability of our filtering method to eliminate the stochastic 
noise present in heated configurations, and make use of the information provided to monitor 
the thermalization process itself. 

In practise, we have concentrated most of our work in pure SU(2) gauge theories. 
Our initial configurations are in all cases Q = 1 instantons of various sizes obtained with 
periodic or twisted boundary conditions. If the size parameter of the instanton is small 
compared to the size of the four-dimensional box no significant change of the instanton 
shape is expected with respect to the infinite volume continuum instanton. 

We have applied heat-bath sweep updates to these smooth configurations. Short 
wavelength stochastic noise comes in very fast and short range observables (as the average 
plaquette value) quickly approach values close to those corresponding to thermalized en- 
sembles. The size of the noise is of order l/y^, but even for relatively large values of /3 
it dominates the value of the action density, and the trace of the initial instanton appears 
to be washed out. This can be clearly observed in figure 9. The left pictures display the 
self-dual part of the action density of a smooth instanton integrated in two directions. The 
center figures display the same quantity after applying Nh heat-bath sweeps to it, at two 
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different values of /3 (/3 = 3 with = 10 for the top figure, and /3 = 8 with = 60 
for the bottom one). Notice the tremendous change in scale which follows from the size 
of stochastic noise, even at these relatively large values of /3. Nevertheless, the original 
instanton is still there. The right figures display the image obtained after applying our fil- 
tering method to the configurations depicted in the center figures. Remarkably, the shape 
looks rather similar to the one of the original instanton, and certainly not noisier. 

Although the results look spectacular we are not satisfied by just obtaining a clean 
smooth picture of the underlying instanton. We aim at having full control of the possible 
distortions introduced in the parameters (position and size) of the initial instanton. This 
will be analysed in the next section, where we will describe the ability of the filtering 
method to extract the stochastic motion performed by these moduli parameters under 
heating. 

On the basis of the perturbative analysis of section 2 we knew that the filtering method 
should act appropriately for configurations which are deformations of solutions of the clas- 
sical equations of motion. However, why should the method still work so well in these 
cases in which the size of the noise is rather large? In the remaining of this section we will 
provide a partial answer to this point. 

For that purpose, we studied the effect of heating on the low-lying spectrum of the 
operator 0^^\ The initial smooth configuration was taken to be an SU(2) instanton of size 
p = 3.5a on a lattice of length Ng = 14, employing both periodic and twisted boundary 
conditions for the reasons explained previously. Starting from these smooth configurations 
we performed Nf^ = 10 heat-bath sweeps with the Wilson action at various values of /?, 
ranging from /3 = 3 to /3 = 30. For each initial configuration and value of /3 we repeated 
the procedure 10 times with independent random seeds, and averaged the values obtained 
for the four lowest eigenvalues of O^^^ over them. The results are displayed in figure 10. 
The continuous line gives the result of a fit to the two smallest eigenvalues up to /3 > 3.3 
for each boundary condition. For twisted boundary conditions we get: 
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0.190(2)^+2.02(7)^ , (4.1) 



A^ = ^ [0.57(1) ^ + 0.24(1)- + 2.3(3)^ 



'/3 ' ^ 

7r2 . , 1 ,1 



(4.2) 



and for periodic boundary conditions we obtain: 



\pbc 1 



0.191(2)^ + 1.94(5)^ , (4.3) 



Af- = - 0.9(2) + 0.213(9)- + 1.98(25)^ . (4.4) 



1 



Notice that as /3 goes to infinity the results approach those obtained for smooth instantons 
in the previous sections. There is a zero-mode (up to lattice cut-off effects) for both 
boundary conditions, and a gap which is volume dependent and strongly dependent on 
boundary conditions. For large values of /?, the lowest eigenvalue is proportional to 1//3 
as predicted by perturbation theory. The leading correction to the higher eigenvalues 



-23- 




0.25 



CO 



0.15 




0.05 



Figure 10. We display, as a function of the four lowest eigenvalues of the O^^'' operator 
for the ensemble of heated instantons described in the text. Left and right plots correspond to 
twisted and periodic boundary conditions respectively. Each data point represents the average of 
the corresponding eigenvalue over a sample of ten configurations obtained by applying Nh heat- 
bath sweeps to a smooth Q = 1 instanton. The errors are the dispersion over the sample set. The 
continuum lines are fits of the data to a third degree polynomial in - see Eqs. (4.1)-(4.4). 



seems to behave in the same way. Perturbation theory allows that deformations lead to 
a modification of the excited eigenvalues proportional to Nevertheless, since the 

noise introduced by the heat-bath method is stochastic and has both signs, we find it quite 
possible that in this case the correction to the eigenvalue is also of order 1//3 as our data 
seem to indicate. Furthermore, it is remarkable that a fit to the /3 dependence, depicted 
as a continuum line in figure 10, shows that the coefficient of the 1//3 term is very much 
compatible for the first and second eigenvalue and both boundary conditions. The same 
applies for the coefficient of the 1//?'^ term. There is no particular reason why there should 
be no l//3^ term as well, but adding it to the fit gives a coefficient which is compatible 
with zero, and a larger value for the P^r degree of freedom. 

Some final comments are in order. The lowest eigenvalue of the oj^^'' operator depends 
mildly on the value of N^- Typically, it starts to rise for the initial heat-bath steps, but 
quickly saturates to a fixed value. The value also depends slightly on the choice of the 
heat-bath action. Taking the improved e = action instead of the Wilson action only 
affects the coefficients by 10%. 

The most remarkable feature shown by figure 10, is that the gap remains approximately 
constant for all values of /3. This is the reason why our method works well in all that range. 
Had level crossing occurred, our filtered mode would have lost the trace of the underlying 
instanton, and failed dramatically. In the following section we will go one step further and 
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investigate how the heat-bath process and the filtering affect the instanton parameters. 
5 Analysing Monte Carlo dynamics 

In this section we wih present the results of a detailed study of the effect of heat-bath 
updates and subsequent filtering upon the collective parameters describing an instanton: 
its position and size parameter p. The philosophy is similar to that presented in the 
previous section but taking spatial care in minimising systematic errors and increasing the 
statistical significance of the results. 

Our initial smooth instanton configuration was generated by e-cooling [20] with e = 
in order to reduce lattice artifacts. Different instantons sizes have been studied, but most 
of our data has been obtained for p = 3.5a. This is a safe intermediate value, far from 
the p = 2a, in which the instanton soon evaporates and the lattice action is strongly 
dependent of the p value, and the larger values where finite volume starts affecting the 
shape and stability. The lattice size was taken as Ng = 16 and we used twisted boundary 
conditions (m = 0, /e = (1, 1, 1)) for the reasons explained before. 

The Monte Carlo updating procedure consisted on a heat-bath algorithm based on 
the classically improved e = lattice action described previously. Two values of /3, 16 
and 8, were used in our analysis. These relatively large values of /3 ensure that the heat- 
bath updating mostly adds stochastic noise to the configuration, but does not destroy it 
or create new structures. Having two values of /3 allows to quantify how much of these 
fluctuations are gaussian, allowing us to extrapolate the results for lower values of /3. 
Finally, we have studied the evolution of the sample for different number of heating steps 
Nh, by applying the filtering procedure to all the noisy configurations at specific values of 
Nh and extracting from the AFM density the value of the instanton parameters, using the 
methodology described in the appendix. 

How should the distribution of instanton parameters look like for the heated sample? 
It is quite clear that, in addition to generating gaussian fluctuations, the heat-bath dy- 
namics should also generate modifications of the moduli parameters characteristic of the 
underlying smooth configuration. In particular, for instantons, both the location and size 
parameter p should behave stochastically as a result of the updating procedure. Up to 
lattice artifacts the position should describe a random walk with no drift term. On the 
other hand, the quantum fiuctuations should not be fiat in p space, because it is precisely 
the quantum fluctuations introduced by the updating procedure what leads to the size 
dependence of the free energy computed by 't Hooft [2]. Although, the updating procedure 
is discrete, we expect that for sufficiently large values of j3 and of the number of updat- 
ing steps the evolution could be well modelled by a Langevin equation. The equation for 
the displacement parameter being purely a Wiener process, with a probability distribution 
satisfying the heat equation. More precisely, if we started with a gaussian distribution, it 
will remain gaussian at all times, with constant mean value and a dispersion behaving as 



a^{t) = {{xi{t) - xmf) = a\Q) + tK^, 



(5.1) 
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Figure 11. Histogram of position displacement at Nh = 40 and /3 = 16. The solid line is the fit to 
a gaussian. 

The slope parameter measures the velocity of the diffusion process and is a characteristic 
of the updating procedure. On general grounds we know, however, that it is proportional to 
the square of the size of the noise. Thus, for our particular case it should go like oc 1//3. 

Let us see how our data behaves with respect to these expectations. For that purpose, 
starting from our smooth instanton of size p/a = 3.5, we generated 7885 independent 
trajectories by applying a sequence of up to 120 heat-bath steps with j3 = 16. We stored 
the configurations after A'^^ =25, 40, 80 and 120 heat-bath sweeps. To the four sets of 7885 
configurations we applied our filtering procedure. In all cases the filtered image showed a 
clear peak close to the position of the initial instanton similar to that displayed in figure 9. 
For each filtered image we extracted the position and size of the corresponding instanton 
following the procedure specified in the appendix. Finally, we performed an statistical 
analysis of the results. In particular, we computed the displacement of the center of the 
instanton obtained from the AFM relative to the position of the initial instanton. The data 
agrees nicely with the expectations of a gaussian distribution. For example, figure 11 shows 
the statistical distribution for = 40 of the displacement putting the 4 components of 
the displacement in the same plot to increase the statistics. The data are compared with 
a gaussian distribution fitted to the data. The pe^' degree of freedom is 1.15. For the 
remaining values of A'^^ similar results were obtained with x^-pdof of order 1. Finally, the 
mean value and dispersion of the distribution were obtained either directly or by a fit to 
a gaussian distribution, with compatible results. The mean value obtained was always 
smaller than 0.01 lattice spacings with no particular A'^^ dependence. On the contrary, the 
dispersion was observed to grow linearly with the number of heat-bath steps. Figure 12 
shows the values of the dispersion as a function of Nh/ (3. The solid line is the result of a 
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Figure 12. Dispersion in displacement as a function of ^ for /? = 16. The solid line is a linear 
fit. 

fit to a straight line having ^ per degree of freedom smaller than 1. The fit is given by: 



with Ani- = 0.0120(3), and = 0.0286(2). The relation of with the slope parameter 
Kx entering Eq. (5.1) is given by NhB^ = t/3Kx. All results are given in units of the lattice 
spacing. 

To test the /3 dependence we repeated the procedure for [3 = 8, = 20, 40 and 
60 heat-bath sweeps, and 3940 trajectories. The results are qualitatively the same as 
described for /3 = 16. The dispersion is fitted to 0.022(3) + 0.032(2)^. Notice that the 
slope is compatible with the one obtained at /3 = 16, in agreement with the hypothesis of 
gaussian fluctuations. One feature that differs from our Langevin model is the fact that 
the intercept is non-zero. It is not surprising that the continuous Langevin description 
does not apply during the first few steps. Indeed, the agreement is expected only in the 
limit P — > oo with N^/ (3 fixed. This seems to affect mostly the intercept. Apparently the 
intercept behaves as 0.19//3. In any case, this small change is the only one which could 
reasonably be attributed to the distortions associated to the filtering method itself. 

Now we proceed to the analysis of the p evolution. In this case we expect a drift due to 
the p-dependence of the free energy. According to the calculation of Ref. [2] the probability 
distribution of the size parameter should behave as 



((Ax,)') = A, + i?.-^ 



(5.2) 




= P 



7/3 



The first factor comes from the scale invariant distribution in collective coordinates, while 
the second one follows from the exponential of the action with a renormalized coupling 
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Figure 13. Dispersion in p as a function of for /3 = 16 and /3 = 8. The solid lines are linear 
fits. 

constant computed at the scale p. Indeed, the value 22/3 follows from the leading-order 
beta function. The free energy is then given by minus the logarithm of this term. Thus, 
the Fokker-Planck equation which describes the stochastic evolution of the distribution in 
p is given by 

In this case, even if one starts with a gaussian distribution in p at time equal zero, it will 
cease to be gaussian at longer times. Nevertheless, numerical integration of this equation 
reveals very small deviations from gaussianity during the range of times corresponding to 
our heat-bath updates. The drift term shows up in the prediction that the mean value 
{p) becomes an increasing function of time. This function is approximately a straight line 
in the appropriate time range. The dispersion, a'^ip) = {p^) — {p)'^, shows also a (very 
approximate) linear growth with time, with a slope given by Kp. One can circumvent even 
these minor deviations by computing (p^). According to the Fokker-Planck equation one 
can easily deduce that this quantity should grow exactly linearly with time: 

{p'){t) = {p')iO) + K'pt . (5.4) 

The corresponding slope should be higher than Kp by a factor that gives a measure of the 
beta function: 

K 7 

Let us now see how the data behaves with respect to our expectations. In figure 13 
we plot the value of o"^(p) for our /3 = 16 and /3 = 8 data, together with a linear fit to the 
data. The slopes are compatible (0.0168(4) for /3 = 16 and 0.0175(16) for /3 = 8). Again 
the intercept scales as 0.36//3. 
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Figure 14. The expectation value of is displayed as a function of the number of heating steps 
over f3. The upper data points corresponds to /3 = 8, and the lower points to /3 = 16. The two 
continuum lines for each j3 value correspond to different linear fits, with either two free parameters 
or with the slope fixed to the predicted value: B'^ = 10i?p/3. 

The most crucial aspect is the behaviour of {p'^){t) which, according to our expecta- 
tions, should grow linearly with a slope which is 10/3 larger than the one of cr'^ip)- In 
figure 14 a fit to the data is performed showing, for /? = 16, a slope of 0.050(4). This 
gives a ratio of 8.9(8)/3 instead of 10/3. This is less than 2 sigma away from the expected 
value, and gives in any case a value incompatible with the one observed for the dispersion 
(T^(p). We can safely claim that we have obtained numerical evidence for the predicted 
dependence of the free energy with p. The data at /? = 8 show the same trend, although 
this time the slope overshoots and the ratio is rather 13(l)/3. All the parameters of the 
linear fit to our data are collected in table 4. The errors quoted include both statistical 
(using jack-knife) and systematic effects coming from correlation upon fit-parameters. 

We have done other checks at different values of p but with smaller statistics. It is 
impossible in this case to determine the drift in the mean value of p, but given that the 
dispersions are affected by much smaller errors, we can extract meaningful values of Kx 
and Kp from the data. In particular we generated 500 trajectories at p/a = 2.85, 3.25 and 
3.784. In all cases we obtain good linear fits for the dispersion in displacement c7^(Ax) 
and in size cr'^{p)- The parameters of the fits are collected in table 4 and are reasonably 
compatible with each other. Although, the range in p is small and the errors are large, we 
can safely conclude that the slopes are not strongly dependent on the value of p. 

Our data predicts the dependence of the displacement and p as a function of the 
number of heating steps. The results can be extrapolated to lower values of /3. Although, 
as P decreases we expect non-linear dependencies in /? it would be very surprising that the 
results change by orders of magnitude. What are the possible implications of these results? 
Suppose one has a Monte Carlo generated configuration with a particular instanton content. 
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B 


Pn /a 


# Confs 


Nh 


/3^x 


Bx 


l3Ao 


Bn 
P 


16 


3.5 


7885 


25, 40, 80, 120 


0.192(7) 


0.0286(4) 


0.357(15) 


0.0168(4) 


8 


3.5 


3940 


20, 40, 60 


0.17(2) 


0.032(2) 


0.36(5) 


0.0175(16) 


16 


3.25 


500 


5, 10, 20, 50, 100 


0.10(1) 


0.031(1) 


0.20(4) 


0.019(2) 


16 


3.78 


500 


5, 10, 20, 40 


0.14(2) 


0.034(2) 


0.28(6) 


0.023(4) 


16 


2.8 


500 


10, 20, 30, 40 


0.06(1) 


0.0314(5) 


0.17(2) 


0.0189(15) 



Table 4. List of fitted parameters to the dispersion in position displacement ((T^(Aa;) = + ^"^^ ) 

and p values {cP'^p) = Ap + ^^^) ■ The relation with the continuum slope parameters and Kp 
is given by: NhB^ = t/3A'^, and N^Bp = tfiKp. 



Applying heat-bath sweeps to this configuration the positions and sizes of instantons will 
change approximately according to the pattern derived in this paper. Eventually either 
by displacement or by growth the instantons will meet anti-instantons and with a finite 
probability annihilate. The last steps could be very fast, but for large initial separations 
most of the time will be employed in the stochastic motion depicted in this paper. Other 
processes are possible. Equal sign instantons might merge into new structures and even 
percolate as a result of the growth in p. In addition, dislocations can be created which 
might eventually develop into new instantons and change the value of the topological charge 
of the configuration. Finally, one can have the reverse phenomenon, since small instantons 
have a finite probability of decreasing in size and disappearing through the holes of the 
lattice. It is clear that only once these processes have acted sufficiently can we approach a 
thermal distribution of the topological structure. The autocorrelation time for this type of 
processes are much larger than those affecting local quantities such as the mean value of 
the plaquette. Assuming that does not depend on p/a, as our data seems to support, 
we can conclude that, in order to reach changes in displacements and sizes which are given 
by appropriate values I in physical units, we need to perform a number heating steps given 
by: 

This grows very fast for large values of /3 due to the presence of the a^(/3) factor in the 
denominator. This is a typical critical slowing down phenomenon. The factor of 1/-B adds 
to the difficulty, since it is of order 30-60 for displacements and sizes respectively. If we 
input the characteristic values of / and (5 one can provide an estimate of the corresponding 
autocorrelation times. 



6 Conclusions and Outlook 

In this paper we have studied a filtering method proposed in Ref . [36] . We have considered 
a simple and precise definition of the filtered profiles, and implemented it on the lattice 
by means of the overlap Dirac operator in the adjoint representation. We have tested the 
method with various classes of initial configurations. First of all, we analysed solutions 
of the classical equations of motion for which the method should work optimally in the 
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continuum. The goal here was to determine the effects of lattice artifacts and the dis- 
tortions obtained in the extraction of the physical parameters of the configuration. For 
the case of instantons these parameters are the position and the size parameter p. Then, 
we proceeded to study smooth configurations which are not classical solutions. In par- 
ticular, we focused on instanton-anti-instanton pairs at varying distances. For all types 
of smooth configurations the method performs nicely until the configuration is too rough 
(small instantons) or very far from a classical configuration (very nearby instantons and 
anti- instantons). Widely separated structures might also give rise to problems, since the 
gap is reduced. Our goal in this study was not only to see that the method works nicely 
with very small distortions, but also to identify possible limitations and problems that can 
arise when applying the method to a typical configuration appearing in the QCD vacuum. 

Next, we proceeded to add quantum fluctuations by applying a series of heat-bath 
sweeps to an initial smooth configuration. To guarantee that the noise does not destroy 
the initial topological structure we considered a lattice inverse coupling (3 larger than 3 
and a small number of heat-bath steps. The results are quite impressive since the filtered 
image looks very smooth and remarkably close to the initial configuration. The variations 
in the parameters of the filtered image with respect to those of the initial configuration 
cannot be entirely attributed to distortions induced by the filtering process. Part of these 
variations are due to genuine stochastic motions of these parameters as a result of the 
heating process. Given the good quality of our results, we embarked in a high statistics 
analysis of these stochastic motions for the case of an initial single instanton configuration. 
The results obtained are quite precise and follow perfectly the expectations of a Langevin- 
equation modelling of this process. In particular, the time evolution of the dispersions 
of the instanton parameters (position and size) can be used to estimate the characteristic 
autocorrelation times of heat-bath dynamics for producing changes in the properties of 
topological structures present in the configuration. Furthermore, we see positive evidence 
of the quantum free energy dependence on p computed many years ago by 't Hooft in per- 
turbation theory. Thus, our data provides a numerical determination of the first coefficient 
of the (3 function of Yang-Mills theory, which is roughly in agreement with the number 
extracted from perturbation theory. 

Given that the method seems to work very well in all the tests carried out in this 
paper, it might indeed become an effective tool to analyse the topological content of Monte 
Carlo generated configurations at typical lattice spacing scales. Ultimately it could allow 
to test the validity of all the proposed scenarios of the Yang-Mills theory vacuum at finite 
and zero temperature. This is, of course, a very ambitious goal which demands a thorough 
and careful analysis. Our initial tests look rather promising. 

Fortunately, the method has other applications. In particular, along the lines pre- 
sented in the previous section, one can use it to study the dynamics induced by the Monte 
Carlo updating procedure on characteristic topological structures. There are two aspects 
involved: the shape of the generated quantum distribution and the rate at which this distri- 
bution is attained. With respect to the first aspect, the case of the single instanton studied 
in this paper is the simplest. Perturbation theory provided long time ago the structure of 
the induced quantum weights for small instantons. However, large instantons, calorons and 
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fractional instantons are other topological structures which have a richer parameter space 
and for which perturbation theory has only partial answers. For example, for the case of 
calorons both masses and positions of the constituent monopoles determine its parameter 
space. What is the dependence of the quantum weights on them? There have been some 
recent results in perturbation theory [59] addressing these issues. Our method can verify 
these results and explore other instances for which perturbative results are not available. 

The second aspect is also very important from the practical viewpoint. Indeed, different 
updating procedures might lead to different autocorrelation times for topological features 
present in the configurations. The number of sweeps necessary to generate a thermalized 
sample depends crucially on that number. Our method can then be used to compare 
different updating procedures and determine the optimal one. This issue is a rather hot 
topic which has generated some debate in the literature, after very large autocorrelation 
times were observed for overall topological quantities such as total charge [42]- [44]. 

We want to add some final comments about possible improvements of the method on 
the methodological side. Some of them could be associated with the lattice implementation, 
by using improved versions of the Wilson-Dirac operator entering in the construction of 
the overlap operator. Other changes can be associated with the definition of the continuum 
operator itself O^^-*. In particular, one could take instead operators of the form —PDKDP. 
A convenient choice of the kernel K can preserve all the nice properties of O^^^ (positivity, 
reality, unique zero-mode for classical solutions, etc) while allowing to minimize corrections 
or increase the gap in the spectrum. For the time being these modifications have proved 
unnecessary, but this possibility might be worthy of future exploration. 
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A Instanton size determination 

As explained in the paper, it is important to establish the effect of the filtering procedure 
and of the heat-bath algorithm on the distribution of instanton parameters: position and 
size p. It is clear that a discretized instanton action density is not identical to the analytic 
continuum expression. We expect hence that different lattice determinations of instanton 
position and sizes differ by values of order a. The differences should decrease as p/a 
increases. In addition to this, there is also a finite volume correction which modifies the 
shape of the instanton as p/L approaches 1. This correction is dependent on the boundary 
conditions. 
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In what follows we will describe the methods we have used in this paper to determine 
the instanton parameters. Although they are affected by lattice corrections, our concern 
here is to provide a quantitative comparison between two different density configurations 
(zero-mode or action density, filtered or not-filtered, heated or not-heated) employing the 
same algorithm for extracting the instanton parameters in both cases. In this way an 
important part of the lattice corrections disappear when taking the difference. 

We have used two different methods to determine the instanton parameters from the 
lattice densities. The purpose of this appendix is to explain both of them and compare 
their results. 

The first method extracts the instanton parameters by comparing the lattice density 
field with that of a continuum instanton at infinite volume, given by: 

^^''^ = {p^ + lx - xfY ' ^^-^^ 

where p and X denote respectively the instanton size and position. For that purpose, 
one fits the lattice density field in the vicinity of its local maximum to Eq. (A.I). One 
of the nice advantages of this method is that it can be applied to configurations having 
several sufficiently dilute instantons. The procedure can then be applied to all of the local 
maxima (peaks) and determine all the parameters of each instanton. Following Ref. [21], 
the method attempts to take into account the finite volume effects and the interaction with 
other instanton-like objects by using an ellipsoid generalization of the instanton ansatz 
given by: 

g(.) = i^(l + ^ ("^-/^^Y ' . (A.2) 

For an exact instanton K = 487r^/p^, and all the p^ are equal. This formula has in total 9 
parameters which can be extracted by solving for S^{x) on the lattice location of the peak 
and its 8 nearest neighbours. From the values of we extract a size parameter given by 
the geometrical average: = JI^Pa*- Note that the p parameter computed in this way 
is independent of the overall normalization of the density, which is instead encoded in the 
value of K. 

The second method employs the instanton profile in one dimension. This is obtained 
by integrating the action density over the three spatial coordinates giving 

The corresponding lattice quantity P^(no) is obtained by summing the lattice action den- 
sity over three-dimensional volumes. The parameters Xq and p are then determined by 
identifying the value of P^(no) at the maximum and its two neighbouring lattice points 
with: 

(P^(no))-^/^=/?fl+ ^""°"/°^ 



(A.3) 



As in the previous method, the value of K is also a function of p, but it is not taken into 
account for its determination. Thus, this method is also insensitive to the normalization 
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Figure 15. For several Q = 1 instanton configurations and twisted boundary conditions, we 
display the relative difference between the two size parameters described in the appendix. 



of the density. The procedure is repeated for all directions and in this way we determine 
the four coordinates of the center, and four different estimates of the p parameter. From 
them we can extract the mean value and the error in the determination of p. 

It is clear that the second procedure is non-local as opposed to the first one. This 
implies that in the case of noisy configurations the statistical fluctuations are reduced. 
On the contrary it fails when applied to an environment in which many instantons are 
present. That is the reason why we prefer to employ both methods instead of sticking to 
a single procedure. In any case, for instantons that are sufficiently smooth, the error in 
the size determination induced by the discretization remains at the few percent level. This 
is illustrated in figure 15 where we compare the size parameters, extracted from the two 
methods described above, for a set of twisted Q = 1 instanton configurations of varying 
sizes. For instantons larger than 2.5a, the difference between the two p determinations 
decreases as p/a increases and remains in all cases below 5%. 
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